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Abstract 

With  increasing  developments  in  computer  technology  and  available  software, 
simulation  is  becoming  a  widely  used  tool  to  model,  analyze,  and  improve  a  real  world 
system  or  process.  However,  simulation  in  itself  is  not  an  optimization  approach. 
Common  optimization  procedures  require  either  an  explicit  mathematical  formulation  or 
numerous  function  evaluations  at  improving  iterative  points.  Mathematical  formulation 
is  generally  impossible  for  problems  where  simulation  is  relevant,  which  are 
characteristically  the  types  of  problems  that  arise  in  practical  applications.  Further 
complicating  matters  is  the  variability  in  the  simulation  response  which  can  cause 
problems  in  iterative  techniques  using  the  simulation  model  as  a  function  generator. 

The  mixed-variable  generalized  pattern  search  with  ranking  and  selection 
(MGPS-RS)  algorithm  for  stochastic  response  problems  is  applied  to  an  external 
simulation  model,  by  means  of  the  NOMADrn  MATLAB  software  package.  Numerical 
results  are  provided  for  several  configurations  of  a  simulation  model  representing  a 
multi-echelon  repairable  problem  containing  discrete,  continuous,  and  categorical 
variables.  Computational  experience  results  are  presented. 
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OPTIMIZATION  OF  A  MULTI-ECHELON  REPAIR  SYSTEM  VIA 


GENERALIZED  PATTERN  SEARCH  WITH  RANKING  AND  SELECTION: 

A  COMPUTATIONAL  STUDY 

I.  Introduction 

1.1  Background  of  the  Problem 

With  the  increasing  developments  in  computer  technology  and  available  software, 
simulation  is  becoming  a  widely  used  tool  to  model,  analyze,  and  improve  a  real  world 
system  or  process.  A  simulation  model  is  often  developed  because  the  system  under  study 
is  so  complex  that  an  analytical  model  either  is  difficult  to  develop  or  cannot  be  formulated. 
Each  simulation  model  can  be  classified  as  one  of  two  types  depending  on  the  presence  of 
random,  or  stochastic,  elements  in  the  model.  Those  models  without  random  elements  are 
called  deterministic  models.  Each  time  it  is  run  with  a  particular  set  of  input  variables, 
a  deterministic  model  yields  a  unique  output,  or  response,  with  no  variation.  Simulation 
models  that  include  some  randomness,  which  means  the  response  changes  randomly  for 
each  run,  are  called  stochastic  models.  This  research  addresses  stochastic  models  with 
responses  that  have  variation,  or  noise,  as  a  result  of  randomness  for  each  run  of  the  model. 

Often  an  analyst  uses  the  simulation  model  as  an  ad  hoc  means  to  optimize  the 
real  system.  Typically  with  this  approach,  the  analyst  sets  the  input  variables,  runs  the 
simulation  for  one  or  more  replications,  and  evaluates  the  response.  The  analyst  updates 
the  input  variables  and  repeats  the  process,  ultimately  finding  a  “best”  solution.  The 
input  variables  yielding  that  solution  are  then  implemented  in  the  real  system.  However, 
the  evaluation  of  the  response  with  a  given  set  of  input  variables  is  complicated  due  to  the 
stochastic  nature  of  the  simulation.  While  the  objective  may  be  to  optimize  an  overall 
average  response,  the  response  for  a  specific  set  of  input  variables  may  only  be  the  result  of 
a  small  number  of  simulation  runs.  One  solution  to  this  problem  is  to  run  the  simulation 
at  each  set  of  input  variables  for  a  large  number  of  replications.  However,  this  is  usually 
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not  practical,  particularly  if  the  number  of  input  variable  combinations  to  be  run  is  large 
or  the  simulation  takes  a  long  time  to  run.  What  is  needed  is  a  procedure  that  can  handle 
the  noisy  response  of  the  simulation  to  determine  which  input  variables  produce  the  “best” 
response. 

The  noisy  response  may  be  modeled  as  an  unknown  response  function  F(x,u)  which 
depends  upon  an  n-dimensional  vector  of  controllable  design  variables  x  €  Rn,  and  the 
vector  cj,  which  represents  random  effects  inherent  to  the  system.  The  objective  function 
/  of  the  optimization  problem  is  the  expected  performance  of  the  system,  given  by 

/(x)=Bp[F(x,w)]=  f  F(x,u)P(<kS),  (1.1) 

Jn 

where  u  G  (1  can  be  considered  an  element  of  an  underlying  probability  space  (fi,  F,  P ) 
with  sample  space  Q.  sigma-field  F .  and  probability  measure  P  (51).  Because  the  re¬ 
sponses  come  from  a  black-box  simulation  which  cannot  be  represented  analytically,  the 
probability  distribution  that  defines  the  response  F(x,u>)  is  assumed  to  be  unknown  but 
can  be  sampled. 

An  optimal  solution  for  either  a  deterministic  or  stochastic  simulation  model  can  be 
difficult  to  obtain.  Since  /  is  usually  unknown  and  analytical  derivatives  are  unavailable, 
classical  optimization  approaches  generally  do  not  apply.  Also,  simulation  runs,  necessary 
for  the  numerical  evaluation  of  /,  may  be  computationally  expensive.  The  presence  of 
noise  only  complicates  matters  because  /  cannot  be  evaluated  precisely.  Statistical  tests 
to  determine  if  one  x  is  better  than  another,  a  requirement  for  many  search  methods,  may 
require  a  large  number  of  repeated  samples.  Additional  complications  arise  when  x  con¬ 
tains  non-continuous  variables,  either  discrete-numeric  ( e.g .  integer-valued)  or  categorical. 
Categorical  variables  are  those  that  can  only  take  on  values  from  a  predefined  list  that  have 
no  ordinal  relationship  to  one  another.  For  example,  a  company  may  have  different  types 
of  materials  used  in  the  manufacture  of  their  products.  The  variables  may  represent  those 
materials  ( i.e .  1  =  steel,  2  =  aluminum,  etc).  The  class  of  optimization  problems  that 
includes  continuous,  discrete- numeric  and  categorical  variables  is  known  as  mixed  variable 
programming  (MVP)  problems  (10). 
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1.2  Purpose  of  the  Research 

Sriver  (50)  developed  a  general  algorithm  for  optimization  of  a  stochastic  system. 
Using  a  combined  generalized  pattern  search  with  ranking  and  selection  approach,  the  al¬ 
gorithm  handles  the  mixed  variable  case,  where  the  input  vector  contains  continuous  and 
discrete/categorical  variables,  with  bound  and  linear  constraints  to  optimize  a  black-box 
simulation  response.  This  thesis  builds  on  the  work  of  Sriver  by  developing  a  simula¬ 
tion  representing  a  real  world  system  to  apply  the  Mixed-Variable  Generalized  Pattern 
Search  with  Ranking  and  Selection  (MGPS-RS)  algorithm  for  stochastic  response  func¬ 
tions.  Dunlap  (24)  incorporated  the  MGPS-RS  algorithm  into  NOMADnr,  an  optimiza¬ 
tion  tool  written  in  MATLAB1  programming  language  by  Abramson  (1)  and  the  software 
used  to  optimize  the  simulation  model.  The  simulation  model  developed  represents  a 
Multi-Echelon  Repairable  System  that  has  wide  use  in  military  and  industrial  environ¬ 
ments. 

1.3  Overview  of  the  Document 

The  next  chapter  contains  a  review  of  the  literature  for  simulation-based  optimization 
strategies,  including  a  brief  survey  of  pattern  search  and  ranking  and  selection,  followed  by 
a  description  of  real  multi-echelon  repairable  systems  and  related  solution  methodologies. 
Chapter  3  presents  the  application  of  the  algorithm  to  a  simulation  built  to  represent  the 
multi-echelon  repair  problem.  Chapter  4  presents  computational  results.  Finally,  Chapter 
5  offers  some  concluding  remarks  and  recommendations  for  further  advancement  of  the 
optimization  via  simulation  using  the  MGPS-RS  stochastic  algorithm. 


1 

MATLAB  is  a  registered  trademark  with  MathWorks. 
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II.  Literature  Review 


Prior  to  applying  the  algorithm  for  the  optimization  of  the  system,  it  is  important  to 
review  the  literature  for  competing  methods  and  examine  the  procedures  used  in  the  algo¬ 
rithm.  Section  2.1  provides  an  overview  of  techniques  and  methods  for  simulation-based 
optimization.  The  approach  used  in  this  thesis,  due  to  Sriver,  uses  a  generalized  pattern 
search  method  modified  to  handle  the  stochastic  nature  of  optimization  primarily  for  its 
independence  from  the  need  for  gradient  information  and  for  its  convergence  theory.  Sec¬ 
tion  2.2  gives  a  history  and  recent  advancements  of  the  pattern  search  class  of  algorithms. 
Sriver’s  work  focused  on  extending  pattern  search  methods  to  stochastic  problems  with 
noisy  responses.  He  accomplished  this  by  augmenting  pattern  search  with  a  ranking  and 
selection  strategy.  Section  2.3  describes  some  general  ranking  and  selection  approaches. 
The  last  section  discusses  the  multi-echelon  repair  problem  and  some  of  the  recent  solution 
methodologies,  including  approaches  that  use  simulation-based  optimization. 

2.1  Simulation-based  Optimization 

There  are  several  papers  that  discuss  the  foundations,  theoretical  developments,  and 
applications  of  a  variety  of  techniques  for  simulation  optimization  (32,  26,  9,  13,  53).  Rank¬ 
ing  and  selection,  described  more  thoroughly  in  Section  2.3,  is  a  popular  methodology,  but 
it  does  not  handle  a  large  number  of  candidate  solutions  and  is  impractical  in  the  case 
of  continuous  variables.  When  applying  these  procedures  to  problems  with  continuous 
variables,  the  variables  must  be  discretized.  The  intervals  are  often  user-defined  and  can 
combinatorically  explode  the  search  space  when  not  appropriately  specified.  Multiple 
comparison  procedures  run  a  number  of  replications  and  make  conclusions  on  a  perfor¬ 
mance  measure  by  constructing  confidence  intervals  (26).  Likewise,  multiple  comparison 
procedures  work  better  when  the  entire  decision  space  is  completely  discrete.  Random 
search  can  deal  with  a  large  number  of  candidate  solutions,  as  well  as  upper  and  lower 
bounds.  However,  because  previous  information  is  not  used  at  each  iteration,  formal  con¬ 
vergence  proofs  for  random  search  methods  are  rare,  especially  with  continuous  variables 
and  noisy  response  (49). 
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Response  surface  methodology  (RSM)  is  a  class  of  procedures  characterized  by  fitting 
a  series  of  regression  models  to  the  responses  from  a  simulation  evaluated  at  several  specific 
design  points,  then  optimizing  the  resulting  regression  function.  RSM  is  a  popular  method 
because  of  its  use  of  well-known  statistical  properties.  In  application  to  simulation-based 
optimization,  much  of  the  research  in  polynomial  based  RSM  prior  to  1990  is  summa¬ 
rized  in  Jacobson  and  Schruben  (32),  in  which  several  improvements  are  discussed  such  as 
screening  for  variable  reduction,  allowance  for  multiple  objectives,  constraint-handling  via 
the  methods  of  feasible  directions  and  gradient  projection,  variance  reduction  via  common 
and  antithetic  pseudorandom  numbers,  and  the  effects  of  alternative  experimental  designs. 
RSM  does  have  drawbacks,  notably  its  lack  of  convergence  and  inability  to  handle  cat¬ 
egorical  variables  (50).  Gradient-based  methods,  such  as  finite  difference,  perturbation 
analysis,  and  likelihood  ratio,  that  estimate  gradients  of  the  objective  function  are  well- 
known  and  widely  used,  but  are  restricted  to  the  continuous  variable  problem.  Stochastic 
approximation  is  also  a  gradient-based  method  that  recursively  estimates  the  gradient. 
This  method  possesses  some  convergence  theory  and  certain  variants  can  be  quite  effi¬ 
cient,  but  like  the  other  gradient-based  methods,  it  is  geared  towards  continuous  variable 
problems.  The  preceding  techniques  can  be  classified  into  two  types:  those  for  use  with 
discrete  input  parameters  and  those  for  use  with  continuous  variables.  None  of  these  are 
able  to  deal  with  the  mixed  variable  problem. 

Much  of  the  simulation  software  available  today  contains  some  sort  of  optimization 
package,  usually  in  the  form  of  a  heuristic  search  (27).  A  search  heuristic  is  a  method 
used  to  solve  a  problem  essentially  by  trial  and  error.  The  procedure  described  at  the 
beginning  in  this  thesis  is  also  an  example  of  a  heuristic.  While  heuristics  often  have  an 
intuitive  justification  and  can  yield  good  solutions,  they  do  not  necessarily  produce  an 
optimal  solution  (56).  Examples  of  search  heuristics  used  in  simulation  software  include 
evolutionary  algorithms  (genetic  algorithms,  evolutionary  strategies  and  evolutionary  pro¬ 
gramming),  scatter  search,  tabu  search,  and  simulated  annealing.  Their  relative  ease  of 
use  and  generality  (they  can  easily  be  adapted  to  mixed-variable  problems  and  require  only 
black-box  response  samples)  have  made  them  popular  choices  for  simulation-based  opti¬ 
mization.  However,  their  application  to  stochastic  problems  has  been  largely  unmodified 
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Table  2.1  Summary  of  techniques  proposed  for  simulation-based  optimization. 


Continuous 

Di: 

Numeric 

screte 

Categorical 

Convergence 

Ranking  &  Selection 

X 

X 

X 

Multiple  Comparison 

X 

Random  Search 

X 

X 

X 

RSM 

X 

X 

Gradient-Based 

X 

X 

Heuristics 

X 

X 

X 

from  their  original  form,  relying  on  inherent  robustness  to  noise,  rather  than  explicitly  ac¬ 
counting  for  noise  (27).  Boesel  et  al.  (15)  provide  a  good  framework  for  simulation-based 
optimization.  Their  work  is  similar  to  the  algorithms  used  in  this  thesis,  particularly  their 
use  of  ranking  and  selection.  However,  they  employ  the  genetic  algorithm  heuristic,  which 
requires  the  continuous  variables  to  be  discretized.  Moreover,  the  user  must  choose  the 
discrete  interval  for  each  continuous  variable. 

Table  2.1  summarizes  this  section,  listing  the  techniques  with  their  appropriate  use 
and  convergence  theory.  A  field  containing  an  “X”  denotes  the  ability  of  the  particular 
method  to  handle  a  certain  type  of  variable  or  shown  to  have  convergence.  The  small 
“x”  signifies  that  random  search  methods  have  some  convergence  theory,  but  not  when 
continuous  variables  are  part  of  the  decision  space. 


2.2  Pattern  Search 

Pattern  search  methods  are  a  class  of  direct  search  methods  for  nonlinear  optimiza¬ 
tion.  The  term  “direct”  means  the  methods  use  minimal  information  about  the  objective 
function,  making  a  direct  comparison  of  objective  function  values  without  the  need  for  ex¬ 
act  derivatives  or  their  approximation.  Since  an  objective  function  from  a  simulation  may 
not  be  easily  computed,  direct  search  methods  apply  well  to  simulation-based  optimiza¬ 
tion.  Also,  because  of  their  broad  applicability,  allowing  for  mixed  variable  input  vectors, 
they  can  be  applied  to  simulations  with  a  variety  of  discrete  and  continuous  parameters. 
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Pattern  search  methods  can  be  traced  back  to  the  work  of  Box  in  industrial  efficiency 
during  the  1950s  (18)  as  a  way  to  create  a  “layman’s”  method,  not  relying  on  the  experi¬ 
mental  designs  and  regression  of  response  methodology.  Additional  direct  search  methods 
were  developed  in  the  1960s,  but  these  algorithms  were  not  formally  shown  to  have  strong 
convergence  theory  (2).  However,  in  the  late  1990s,  Torczon  introduced  pattern  search 
as  a  generalization  of  several  existing  methods  and  established  convergence  theory  for  the 
entire  class  of  algorithms  (54).  Torczon’s  paper  was  significant  in  that  it  established  a 
global  first-order  convergence  theory  without  ever  explicitly  computing  or  approximating 
derivatives. 

2.2.1  Generalized  Pattern  Search.  Pattern  search  algorithms  are  defined  through 
a  finite  set  of  directions  used  at  each  iteration.  The  direction  set  and  a  step  length  pa¬ 
rameter  are  used  to  construct  a  conceptual  mesh  centered  about  the  current  iterate  (the 
incumbent).  Trial  points  are  selected  from  this  discrete  mesh,  evaluated,  and  compared 
to  the  incumbent  in  order  to  select  the  next  iterate.  If  an  improvement  is  found  among 
the  trial  points,  the  iteration  is  declared  successful  and  the  mesh  is  retained  or  coarsened; 
otherwise,  the  mesh  is  refined  and  a  new  set  of  trial  points  is  constructed.  Torczon  proved 
that,  for  a  continuously  differentiable  function  /,  a  subsequence  of  the  iterates  {£&}  pro¬ 
duced  by  the  generalized  class  of  methods  converges  to  a  stationary  point  of  /  by  showing 
that  the  mesh  size  (step  length)  parameter  becomes  arbitrarily  small. 

The  mesh  is  defined  by  a  finite  set  of  directions  that  must  be  sufficiently  rich  to  ensure 
that  at  least  one  of  them  is  a  direction  of  descent,  provided  that  the  current  iterate  is  not  a 
stationary  point.  Lewis  and  Torczon  (36)  applied  the  theory  of  positive  linear  dependence 
(22)  to  establish  criteria  for  such  a  set  of  directions.  Specifically,  the  set  of  directions  form 
a  positive  spanning  set  for  Rn,  which  is  defined  as  a  set  for  which  any  vector  in  Rn  can 
be  expressed  as  a  nonnegative  linear  combination  of  these  directions.  Typically  this  set 
forms  a  positive  basis,  which  is  the  smallest  proper  subset  of  a  positive  spanning  set  that 
still  positively  spans  W1.  A  positive  basis  contains  between  n+ 1  (a  minimal  set)  and  2 n  (a 
maximal  set)  elements  (22).  Therefore,  the  worst  case  number  of  trial  points  per  iteration 
can  be  bounded  to  n  +  1  points  by  an  appropriately  constructed  direction  set. 
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2.2.2  Bound  and  Linear  Constraints.  Lewis  and  Torczon  extend  the  results  of 
(54)  to  problems  with  bound  constraints  (37)  and  a  finite  number  of  linear  constraints  (38). 
In  these  situations,  the  set  of  positive  spanning  directions  used  in  the  algorithm  must  be 
chosen  such  that  they  conform  to  the  geometry  of  the  nearby  constraint  boundaries.  Using 
this  construct,  at  least  one  direction  in  the  positive  spanning  set  must  be  a  feasible  descent 
direction,  unless  the  current  iterate  is  already  a  stationary  point. 

Audet  and  Dennis  (11)  devise  an  alternative  version  of  pattern  search  for  bound 
and  linearly  constrained  problems,  along  with  a  new  convergence  theory  based  on  the 
nonsmooth  calculus  of  Clarke  (20)  that  generalizes  previous  results.  Second-order  behavior 
is  described  in  (3).  Audet  and  Dennis  explicitly  separate  the  evaluation  of  points  into  two 
distinct  steps,  an  optional  search  and  a  local  POLL  step.  The  step  allows  the  user  to  define 
a  search  strategy  to  seek  an  improved  mesh  point.  The  SEARCH  step  contributes  nothing 
to  the  convergence  theory,  but  allows  the  user  to  apply  any  finite  heuristic  to  increase  the 
efficiency  and  possibly  affect  the  quality  if  a  correct  search  is  chosen  (17).  Examples  for 
the  use  of  this  approach  include  randomly  selecting  a  space-filling  set  of  points  using  Latin 
hypercube  design  or  applying  a  few  iterations  of  a  genetic  algorithm  (50).  For  problems 
with  computationally  expensive  functions,  the  SEARCH  step  is  often  used  to  construct 
inexpensive  surrogate  functions  and  then  optimize  the  resulting  surrogate  problem  ( e.g ., 
see  (17)).  Dunlap  (24)  studied  the  use  of  surrogates  in  pattern  search  methods,  applied 
to  mixed  variable  stochastic  problems.  The  poll  step  evaluates  specific  points  on  the 
mesh,  referred  to  as  the  poll  set,  that  are  adjacent  to  the  current  iterate  with  respect  to 
the  current  set  of  positive  spanning  directions. 

2.2.3  Nonlinear  Constraints.  Audet  and  Dennis  (11)  extend  their  approach  to 
nonlinear  constraints  by  implementing  a  filter  method  (25),  which  accepts  new  iterates  if 
the  usual  improvement  in  objective  function  is  found,  but  also  if  an  aggregate  constraint 
violation  function  is  reduced.  Lewis  and  Torczon  (39)  give  an  alternate  approach  where 
a  bounded  subproblem,  formed  from  an  augmented  Lagrangian  function  (21),  is  solved 
approximately  using  pattern  search.  Motivated  by  weaknesses  in  the  convergence  theory 
of  the  filter  pattern  search  method  (11),  Audet  and  Dennis  (12)  introduced  a  new  class 
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of  algorithms,  called  Mesh  Adaptive  Direct  Search  (MADS),  which  generalizes  GPS  by 
allowing  more  flexibility  in  the  selection  of  directions.  In  fact,  MADS  has  been  proved 
to  converge  to  second-order  stationary  points,  and  even  local  minimizers  of  general  set- 
constrained  nonlinear  optimization  problems  (5).  MADS  uses  a  barrier  method,  replacing 
the  filter  method,  that  assigns  a  value  of  +oo  to  infeasible  iterates  without  ever  evaluating 
their  objective  function. 

2.2.4  Mixed-Variable  Generalized  Pattern  Search.  Audet  and  Dennis  (10)  pro¬ 
vided  a  framework  for  mixed-variable  problems  with  bound  and  linear  constraints  by  in¬ 
cluding  discrete  neighborhood  sets  in  the  definition  of  the  mesh.  Their  algorithm  was 
applied  successfully  in  (35)  to  the  design  of  a  thermal  insulation  system.  In  the  mixed 
variable  case,  the  POLL  step  is  augmented  with  a  search  of  points  in  a  user-defined  set 
of  discrete  neighbors.  If  the  poll  step  is  unsuccessful  in  finding  an  improved  solution, 
an  extended  poll  step  is  performed,  in  which  a  poll  is  performed  around  any  discrete 
neighbor  that  has  an  objective  function  value  sufficiently  close  to  that  of  the  incumbent 
( i.e .  within  a  tolerance  £  >  0).  This  algorithm  allows  extension  of  the  convergence  theory 
to  the  mixed  variable  domain  but  incurs  a  cost  of  more  function  evaluations,  particularly 
if  the  user  allows  a  large  number  of  discrete  neighbors.  Abramson  (2)  extended  the  work 
done  by  Audet  and  Dennis  by  allowing  nonlinear  constraints  in  the  mixed-variable  case 
through  the  use  of  a  filter.  This  work  was  applied  successfully  in  (4)  to  the  design  of  a 
load-bearing  thermal  insulation  system. 

The  following  description  by  Sriver  (50)  explains  the  definition  of  the  mesh  and  poll 
set  developed  by  Audet  and  Dennis  (10): 

“A  set  of  positive  spanning  directions  Dl  is  constructed  for  each  unique 
combination  *  =  1,2,...,  *max,  of  values  that  the  discrete  variables  may  take, 
i.e., 

Dl  =  GiZi,  (2.1) 

where  Gi  6  RnCxn<  is  a  nonsingular  generating  matrix  and  Z\  £  Zn  *\D  I .  The 
mild  restrictions  imposed  by  (2.1)  are  necessary  for  the  convergence  theory. 

The  mesh  is  then  formed  as  the  direct  product  of  Qd  with  the  union  of  a  finite 
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number  of  meshes  in  0C,  i.e. 


'-max  r  i  •  i  a 

Mk{xk )  =  Qdx  (J  \  x%  +  A  kDlz  eOc:z£Zl+  1  L  .  (2.2) 

i= 1  ^  > 

At  iteration  k,  let  D),  C  Dl  denote  the  set  of  poll  directions  corresponding 
to  the  ith  set  of  discrete  variable  values  and  define  Dk  =  D\. .  The  poll  set 

is  defined  with  respect  to  the  continuous  variables  centered  at  the  incumbent 
while  holding  the  discrete  variables  constant.  Its  form  is 

Pk(xk )  =  {xk  +  A k(d,  0)  £  0  :  d  £  Dk}  (2.3) 

for  some  1  <  i  <  imax,  where  ( d ,  0)  denotes  the  partitioning  into  continuous 
and  discrete  variables;  0  means  the  discrete  variables  remain  unchanged,  i.e., 
xk  +  A k(d,  0)  =  (xck  +  A kd,  xf).” 

Within  the  GPS  framework,  mixed  variables  are  incorporated  through  the  use  of 
discrete  neighbors  A f(xk)  at  each  point  xk  in  the  domain.  The  points  in  N(xk)  include  the 
current  point  xk  and  other  points  that  differ  in  at  least  one  of  the  discrete  variables.  For 
example,  if  the  discrete  variables  are  defined  as  integers,  a  neighborhood  structure  may  be 
defined  by  holding  the  continuous  variables  constant  and  changing  only  one  of  the  discrete 
variables  by  a  single  unit,  i.e., 

N{xk)  =  {y  £  0  :  yc  =  x%,  yd  -  xf  <  1}.  (2.4) 

However,  if  the  discrete  variables  are  categorical,  then  this  neighbor  set  may  not  be  well- 
defined.  For  example,  in  a  manufacturing  process,  a  categorical  variable  might  be  material 
type,  in  which  case,  the  norm  function  is  not  well-defined,  since  there  is  no  measure  of 
distance  for  this  nonnumeric  variable.  In  this  case,  changing  the  material  type  from  one 
designated  by  a  “1”  to  one  designated  by  a  “3”  may  just  as  valid  a  changing  it  to  “2” .  Thus 
for  categorical  variables,  the  discrete  neighbor  set  must  be  defined  by  the  user.  It  should 
also  be  noted  that  a  change  in  a  discrete  variable  may  force  an  accompanying  change  to 
the  continuous  variables.  As  an  example,  if  a  continuous  variable,  such  as  thickness,  were 
associated  with  each  material,  then  a  discrete  neighbor  would  be  of  a  different  material 
type  than  the  current  point,  but  it  might  also  have  a  different  thickness. 
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2.2.5  Generalized  Pattern  Search  with  Noisy  Response.  The  use  of  generalized 
pattern  search  applied  to  stochastic  optimization  is  limited.  Trosset  (55)  analyzed  con¬ 
vergence  in  the  unconstrained,  continuous  case  by  viewing  the  iterates  as  a  sequence  of 
binary  ordering  decisions.  For  A&  =  f{X^)  —  f(Y),  where  Xk  is  the  current  iterate  and  Y 
is  a  trial  point  from  the  mesh,  the  statistical  hypothesis  test, 


H0  :  Afc  <  0  (2.5) 

Hi  :  Afc  >  0, 

is  conducted,  in  which  Y  is  accepted  as  the  new  iterate  if  the  null  hypothesis  Ho  is  rejected. 

The  test  is  subject  to  Type  I  and  Type  II  errors.  A  Type  I  error  is  made  if  Ho  is  rejected 

when  it  is  actually  true  and  occurs  with  probability  a,  the  significance  level  of  the  test. 
A  Type  I  error  would  select  a  new  iterate  incorrectly.  A  Type  II  error  is  made  if  Hq 
is  accepted  when  H\  is  true  and  occurs  with  probability  (3.  A  Type  II  error  would  not 
update  the  iterate  with  the  new,  better  point.  The  number  of  Type  I  errors  can  be 
controlled,  ensured  (with  probability  equal  to  one)  to  be  a  finite  number,  by  selecting  a 

OO 

sequence  of  significance  levels  {«&}  such  that  ak  <  oo.  In  addition,  let  {A&}  be  a 

k= 1 

sequence  of  alternatives  satisfying  A j.  >  0,  A*,  =  o( A*.),  and  A*,  — >  0  that  require  power 

OO 

1  —  (3k  when  conducting  the  test  in  (2.5).  Choosing  a  sequence  {/3fc}  such  that  Pk  <  00 

k=  l 

ensures  a  finite  number  of  Type  II  errors  when  A&  >  A&.  Trosset  claims  that  a  sequence 
of  iterates  from  a  GPS  algorithm  can  be  shown  to  converge  almost  surely  to  a  stationary 
point  of  /  but,  in  practice,  would  require  a  very  large  number  of  samples  to  guarantee 
convergence  (55).  He  shows  through  a  power  analysis  that  the  number  of  samples  per 
iteration  grows  faster  than  the  squared  reciprocal  of  the  mesh  size  parameter.  A  power 
analysis  is  a  statistical  technique  to  determine  the  required  sample  size  to  guarantee  a 
probability  (1  —  /?)  of  rejecting  Ho  when  H\  is  true.  Also,  Ouali  et  al.  (43)  applied 
multiple  repetitions  of  generalized  pattern  search  directly  to  a  stochastic  simulation  model 
to  seek  minimum  cost  maintenance  policies  where  costs  were  estimated  by  the  model. 
Sriver  (50)  was  able  to  overcome  the  problem  highlighted  by  Trosset  with  the  use  of  an 
indifference  zone  in  ranking  and  selection  procedures. 
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2.3  Ranking  and  Selection 

Ranking  and  selection  (R&S)  procedures  are  “statistical  methods  specifically  devel¬ 
oped  to  select  the  best  system,  or  a  subset  of  systems  that  includes  the  best  system,  from 
a  collection  of  competing  alternatives”  (28).  The  following  overview  of  R&S  procedures 
details  their  use  in  iterative  search  routines  applied  to  stochastic  optimization  via  simula¬ 
tion. 

Indifference  zone  and  subset  selection  are  the  two  main  topics  in  R&S  methods  (26). 
Indifference-zone  procedures  guarantee  a  solution  within  5  >  0  above  the  true  best  solution 
with  user-specified  probability  (1  —  a),  where  a  €  [0, 1].  The  parameter  <5,  which  represents 
a  measure  of  tolerance  known  as  the  indifference  zone,  is  called  the  indifference  zone 
parameter.  R&S  procedures  collect  response  samples  from  the  alternatives  using  a  single 
stage  or  multiple  stages  of  sampling,  check  a  certain  stopping  criteria,  then  either  continue 
sampling  or  stop  and  select  the  alternative  with  the  smallest  response  estimate  in  the  final 
stage  (49).  The  original  procedure  by  Bechhofer  (14)  determines  the  number  of  samples 
required  of  each  iteration  beforehand,  or  a  priori,  according  to  a  tabular  value  related 
to  the  user-defined  values  for  5  and  a.  A  potential  drawback,  especially  in  the  area  of 
simulation,  is  that  Bechhofer’s  method  assumes  that  the  variance  in  the  response  samples 
is  known  and  equal  across  all  alternatives. 

Dudewicz  and  Dalai  (23)  and  Rinott  (46)  extended  the  approach  to  apply  to  problems 
with  unknown  and  unequal  variances  in  response  samples.  They  used  a  two-stage  process, 
in  which  an  initial  stage  of  sampling  is  done  to  estimate  the  variances,  which  are  then  used 
to  determine  the  number  of  samples  needed  in  the  second  stage  to  ensure  the  probability 
of  correct  selection.  Subset  selection  aims  to  select  a  subset  of  at  most  m  points  and 
guarantees  that  this  subset  contains  at  least  the  best  solution,  with  probability  at  least 
(1  —  a)  (41).  Subset  selection  is  more  useful  when  the  number  of  candidates  is  large. 

To  define  the  requirements  for  a  general  indifference-zone  R&S  procedure,  Sriver  (50) 
gives  the  following  formulation 


“...consider  a  finite  set  {Xi,  X2,  ■  ■  ■ ,  X,nc}  of  nc  >  2  candidate  design 
points.  For  each  i  =  1,2 ,...,nc,  let  /)  =  f(Xf)  =  E[F(Xi,u)\  denote  the 


2-9 


true  objective  function  value.  The  fi  values  can  be  ordered  from  minimum  to 
maximum  as, 

f  [1]  —  f  [2]  —  '  '  '  —  f[nc]'  (2-6) 

The  notation  1m  indicates  the  candidate  with  the  ith  best  (lowest)  true  ob¬ 
jective  function  value.  If  at  least  one  candidate  has  a  true  mean  within  5  of 
the  true  best,  i.e.  /u  —  /m  <  5  for  some  5  >  0  and  i  >  2,  then  the  procedure 
is  indifferent  in  choosing  Xni  or  lu  as  the  best.  The  probability  of  correct 
selection  (CS)  is  defined  as 

P{CS}  =  P  {select  X{1]  \  ,/w  -  f{1]  >  5,  i  =  2, . . . ,  nc]  >  1  -  a,  (2.7) 

where  5  >  0  and  a  £  (0, 1)  are  user  specified.  A  random  selection  of  the 
candidates  guarantees  at  least  P{CS}  =  so  the  significance  level  must 
satisfy  0  <  a  <  1  —  rV” 

When  a  candidate  number  is  large,  traditional  multi-stage  indifference-zone  proce¬ 
dures  may  prescribe  too  many  simulation  runs  because  they  are  based  on  the  least  favorable 
configuration  assumption.  That  is,  the  best  candidate  has  a  true  mean  exactly  5  better 
than  all  other  candidates,  which  are  all  tied  for  the  second  best  (53).  As  a  result,  the  pro¬ 
cedures  can  call  for  an  unnecessarily  high  number  of  samples  in  the  final  stage  to  guarantee 
that  (2.7)  holds. 

Two  recent  directions  in  R&S  research  reflect  attempts  to  address  this  issue.  The 
first  approach  has  been  to  combine  R&S  with  some  type  of  search  strategy  to  explore 
a  large  solution  space.  Olafsson  (42)  and  Pichitlanrken  and  Nelson  (45)  each  introduce 
an  iterative  technique  that  combines  R&S  with  an  adaptive  sampling  algorithm  known 
as  nested  partitioning  (NP),  which  is  used  to  search  the  feasible  space  of  (possibly  large) 
combinatorial  problems  for  a  global  optimum.  The  techniques  use  discrete  time  Markov 
chains  to  show  almost  sure  convergence  to  a  global  optimum  of  the  discrete  variable  space. 

Heuristic  search  methods  with  incorporated  R&S  procedures  have  recently  been  de¬ 
veloped  for  stochastic  optimization.  Ahmed  and  Alkhanris  (8)  describe  and  analyze  a 
globally  convergent  algorithm  for  optimization  over  a  discrete  domain  that  uses  R&S  pro¬ 
cedures  within  simulated  annealing  .  Boesel  et  al.  (16)  and  Hedlund  and  Mollaghasenri 
(30)  apply  R&S  procedures  to  genetic  algorithms. 

The  second  approach  for  resolving  the  high  sample  number  required  for  a  large  set 
of  candidates  combines  the  two  topics  of  subset  selection  and  indifference-zone  selection. 
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These  procedures  identify  and  eliminate  inferior  solutions  and  then  select  the  best  from 
the  remaining  candidates.  Nelson  et  al.  (41)  present  a  general  theory  and  procedure  that 
balances  computational  and  statistical  efficiency.  This  approach  maintains  a  probability 
guarantee  for  selecting  the  best  solution  when  using  the  combined  technique.  Kim  and 
Nelson  (33)  and  Goldsman  et  al.  (28)  present  efficient  fully  sequential  indifference-zone 
techniques  that  eliminate  alternatives  deemed  inferior  as  sampling  progresses. 

Categorical  and  discrete  variables  are  readily  handled  by  modern  R&S  techniques 
since  all  design  alternatives  are  determined  a  priori  and  corresponding  variable  values 
can  be  set  accordingly.  However,  R&S  procedures  have  difficulty  with  a  large  number  of 
solutions.  The  existing  provably  convergent  techniques  (8,  42,  45)  that  combine  R&S  with 
adaptive  search  currently  address  entirely  discrete  domains.  Continuous  variables  can  be 
dealt  with  through  discretization,  but  depending  on  the  interval  chosen,  this  can  cause  a 
combinatorial  explosion  of  the  search  space  and  an  increase  in  computational  expense. 

To  improve  on  the  implementation  by  Trosset  and  applicability  to  the  general  case, 
Sriver  uses  a  ranking  and  selection  procedure  to  identify  a  new  iterate.  Sriver  (50)  lists 
the  specific  advantages  to  include: 

•  It  is  amenable  to  parallelization  techniques  since  several  trial  solutions  can  be  con¬ 
sidered  simultaneously  in  the  selection  process  rather  than  only  two  (incumbent  and 
candidate) . 

•  R&S  procedures  detect  the  relative  order,  rather  than  generate  precise  estimates,  of 
the  candidate  solutions.  This  is  generally  easier  to  do  (27)  and  provides  computa¬ 
tional  advantages. 

•  Selection  error  is  limited  to  Type  II  error  only,  i.e.,  making  an  incorrect  selection 
of  the  best  candidate;  Type  I  error  is  eliminated  based  on  the  assumption  of  a  best 
system  among  the  candidates. 

•  The  use  of  an  indifference  zone  parameter  can  be  easily  and  efficiently  adapted  for 
algorithm  termination. 

Three  such  procedures  were  selected  for  use  in  MGPS-RS:  Rinott’s  two-stage  proce¬ 
dure  (46),  a  screen- and-select  (SAS)  procedure  of  Nelson  et  al.  (41),  and  the  Sequential 
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Selection  with  Memory  (SSM)  procedure  of  Pichitlamken  and  Nelson  (44).  The  SSM 
procedure  is  implemented  by  Dunlap  (24)  in  the  NOMADrn  software  (1)  and  explained 
further  in  Chapter  3. 

2-4  Multi-Echelon  Repairable  Systems 

2-4-1  Basic  Problem.  Multi-echelon,  or  multilevel,  problems  exist  in  many  are¬ 
nas.  One  such  area  that  has  considerable  interest  is  in  the  design  and  performance  of 
maintenance  systems  for  a  repairable  item.  The  general  problem  to  be  investigated  is  the 
determination  of  the  optimal  spare  levels  and  repair  channels  in  a  maintenance  system,  in 
which  a  finite  number  of  items  is  desired  to  be  operational  at  any  given  time,  and  in  which 
queuing  can  occur  at  the  repair  facilities  if  all  channels  are  busy  (6).  The  system  usually 
consists  of  one  or  more  bases  at  the  lowest  level  (or  first  echelon),  one  or  more  depots 
at  the  highest  level,  and  any  number  (to  include  zero)  of  intermediate  levels  in  between. 
Multi-echelon  systems  can  take  on  many  forms  to  include  number  of  levels,  number  of  fa¬ 
cilities  at  each  level,  number  of  machines  in  the  system  that  depend  on  the  use  of  the  item, 
scheduling,  resupply  strategies  (heirachal  and/or  lateral),  transportation  delays,  and  can¬ 
nibalization  or  condemnation.  It  is  because  of  this  complexity  that  multi-echelon  models 
are  often  analyzed  through  simulation  (34). 

2-4-2  Optimization  of  Multi-Echelon  Systems.  Optimization  of  multi-echelon 
models  have  been  approached  both  analytically  and  through  simulation.  Sherbrooke 
developed  the  METRIC  model  to  minimize  expected  backorders,  which  is  equivalent  to 
maximizing  availability  when  no  cannibalization  of  parts  is  assumed  (47).  Using  the  as¬ 
sumption  of  an  infinite  number  of  repair  channels,  Sherbrooke  made  use  of  Palm’s  Theorem 
to  calculate  steady-state  probability  distributions  for  the  number  of  units  due  in  from  re¬ 
pair.  Gross  et  al.  (29)  relaxed  that  assumption  and,  using  Markovian  properties  of  the 
exponential  distribution,  formulated  the  expression  for  machine  availability  in  terms  of  the 
decision  variables,  number  of  spares  and  number  of  repair  channels  at  each  facility.  More 
recent  applications  in  this  area  have  been  based  on  simulation-based  optimization.  Chris- 
sis  and  Gecan  (19)  use  a  direct  search  technique  to  iteratively  seek  an  optimal  solution  to 
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the  multi-echelon  system.  Ahmed  et  al.  (7)  present  an  integrated  approach  of  simulated 
annealing  and  simulation  to  determine  the  design  parameters  of  a  multi-echelon  repairable 
item  inventory  system.  Kochel  and  Nielander  (34)  use  a  genetic  algorithm  to  optimize  a 
five-level  simulation  model.  Alkhamis  and  Ahmed  (6)  use  the  results  from  the  analytical 
solution  by  Gross  et  al.  (29)  to  validate  their  evolutionary  technique  of  particle  swarm 
optimization.  Those  results  are  also  used  in  the  evaluation  of  the  procedures  used  in  this 
thesis. 

2.5  Summary 

This  chapter  has  provided  an  overview  of  the  literature  pertaining  to  this  study.  Sec¬ 
tion  2.1  reviewed  the  various  approaches  to  simulation-based  optimization.  Each  method 
was  reviewed  to  highlight  the  ability  to  handle  different  types  of  variables  and  convergence 
properties.  The  ideas  of  pattern  search  and  ranking  and  selection  used  in  the  MGPS-RS 
algorithm  were  outlined  in  Sections  2.2  and  2.3,  respectively.  Finally,  the  multi-echelon 
repair  problem  was  introduced  in  Section  2.4,  and  some  previous  attempts  at  optimizing 
this  type  of  problem  were  discussed.  The  next  chapter  further  details  the  elements  of  the 
MGPS-RS  algorithm  and  its  use  in  the  NOMADrn  software.  It  also  explains  the  specifics 
of  the  multi-echelon  system  models  used  in  the  application  of  the  MGPS-RS  algorithm. 
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III.  Methodology 

This  chapter  presents  the  methodology  used  in  the  optimization  of  the  multi-echelon  repair 
system.  The  first  section  explains  the  mixed  variable  generalized  pattern  search  with 
ranking  and  selection  algorithm  used  to  optimize  the  system.  Section  3.2  provides  details 
of  the  model  formulation  for  the  multi-echelon  repair  problem.  The  last  section  discusses 
the  simulation  model  and  integration  with  the  optimization  code. 

3.1  MGPS-RS 

When  using  simulation  as  a  black-box  generator  for  objective  function  values,  ob¬ 
viously  the  true  values  in  (2.6)  are  not  available.  Therefore,  it  becomes  necessary  to 
use  samples  of  the  simulation  response  F  to  create  estimates.  Let  nc  be  the  number  of 
candidates  in  the  candidate  set  C.  For  each  i  =  1,2, .. .  ,nc,  let  st  be  the  total  number 
of  replications  and  let  =  {F^Y),,)}*!^  be  the  set  of  responses  obtained  through 

simulation,  where  is  the  input  vector  for  design  i  and  replication  s.  Then  for  each 
i  =  1,2 , ...  ,nc,  the  sample  mean  F)  is  computed  as 

Is' 

Fi  =  ~Y,Fis  l3-1) 

8  =  1 

and  is  an  estimator  for  f% .  These  sample  means  may  be  ordered  in  the  same  manner  as  the 
true  responses  in  (2.6).  The  ranking  and  selection  procedure  determines  the  ordering  of  the 
candidates  with  the  ith  best  estimated  response  value  denoted  by  Yu  G  C.  The  candidate 
with  the  lowest  mean  response,  Yrn ,  is  chosen  as  the  best  point.  A  general  R&S  algorithm, 
as  given  by  Sriver  (50),  is  shown  in  Figure  3.1,  where  the  input  parameters  include  a 
candidate  set  C,  significance  level  a ,  and  indifference  zone  5.  Procedure  RS((7,  a,  S) 
returns  the  best  candidate  Ym. 
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Procedure  RS(C,  a,  5) 

Inputs:  A  set  C  =  {Pi,  Y2, _ ,  Ync}  of  candidate  solutions,  significance  level  a,  and 

indifference  zone  parameter  5. 

Step  1 :  For  each  candidate  Yq.  use  an  appropriate  technique  to  determine  the  number  of 
samples  s*  required  to  meet  the  probability  of  correct  selection  guarantee,  as  a  function 
of  a,  6  and  response  variation  of  Yq. 

Step  2:  Obtain  sampled  responses  Fjs ,  i  =  l,...,nc  and  s  =  1 ,  ...,Sj.  Calculate  the 
sample  means  F)  based  on  the  s*  replications  according  to  (3.1).  Select  the  candidate 
associated  with  the  smallest  estimated  sample  mean,  Fni  as  having  the  h-near-best  mean. 

Return:  Ppj 

Figure  3.1  MGPS-RS  Algorithm  for  Stochastic  Optimization 

Sriver  incorporates  this  R&S  procedure  into  the  MGPS  algorithm  for  deterministic 
mixed  variable  optimization,  due  to  Audet  and  Dennis  (10).  The  MGPS-RS  algorithm 
of  Sriver  (50)  is  presented  in  Figure  3.2.  The  binary  comparison  of  the  incumbent  and 
trial  points  used  in  the  deterministic  case  is  replaced  with  Procedure  RS((7,  a,  8).  The 
algorithm  can  use  any  specific  R&S  procedure,  as  long  as  it  satisfies  the  probability  of 
correct  selection  guarantee  given  in  (2.7). 
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MGPS-RS  Algorithm  for  Stochastic  Responses 

Initialization:  Set  the  iteration  counter  k  to  0.  Set  the  R&S  counter  r  to  0.  Choose  a 
feasible  starting  point,  Xq  £  0.  Set  Ao  >  0,  £  >  0,  ao  £  (0, 1),  and  S o  >  0. 

1.  Search  step  (optional):  Employ  a  finite  strategy  to  select  a  subset  of  candidate 
solutions,  Sk  C  Mk(Xk)  defined  in  (2.2)  for  evaluation.  Use  Procedure  RS(SfcU{Ah.}, 
ar,  Sr)  to  return  the  estimated  best  solution  Ym  £  St  U  {Xk}-  Update  ar+i  <  cxr, 
<5r+i  <  Sr,  and  r  =  r  +  1.  If  Ym  A  Xk,  the  step  is  successful,  update  Xk+i  =  Ym, 
Afc+i  >  Ak*  and  k  =  k  +  1  and  repeat  Step  1.  Otherwise,  proceed  to  Step  2. 

2.  Poll  step:  Set  extended  poll  trigger  £k  >  £.  Use  Procedure  RS (Pk{Xk)  U  Af(Ah-), 
ar,  5r)  where  Pk(Xk)  is  defined  in  (2.3)  to  return  the  estimated  best  solution  Ym  £ 
Rfc(Afc)  U N(Xk).  Update  ar+ \  <  ar,  5r_|_i  <  5r,  and  r  =  r+  1.  If  Ym  /  Xk,  the  step 
is  successful,  update  Xk+\  =  Ym,  A^+i  >  Ak*  and  k  =  k  +  1  and  return  to  Step  1. 
Otherwise,  proceed  to  Step  3. 

3.  Extended  poll  step:  For  each  discrete  neighbor  Y  £  M(Xk)  that  satisfies  the 
extended  poll  trigger  condition  F(Y)  <  F(Xk)  +  £k,  set  j  =  1  and  Yk  =  Y  and  do 
the  following. 

(a)  Use  Procedure  RS(Pfc(Y^),  ar,  5r)  to  return  the  estimated  best  solution  Ym  £ 
Pk(Yk).  Update  ar+\  <  ar,  Sr+ 1  <  Sr,  and  r  =  r  +  1.  If  Ym  /  Yk,  set 
Y^+1  =  Ym  and  j  =  j  +  1  and  repeat  Step  3a.  Otherwise,  set  Zk  =  Yj!  and 
proceed  to  Step  3b. 

(b)  Use  Procedure  RS(A*.  U  Zk,  ari  5r)  to  return  the  estimated  best  solution  Ym  = 
Xk  or  Ym  =  Zk.  Update  ar+\  <  ar,  dr+i  <  dr,  and  r  =  r  +  1.  If  Ym  =  Zk,  the 
step  is  successful,  update  Xk+\  =  Ym,  A^+i  >  Ak*  and  k  =  k  +  1  and  return 
to  Step  1.  Otherwise,  repeat  Step  3  for  another  discrete  neighbor  that  satisfies 
the  extended  poll  trigger  condition.  If  no  such  discrete  neighbors  remain,  set 
Xk+i  =  Xk,  Afc+1  <  Afc*  and  k  =  k  +  1  and  return  to  Step  1. 

*NOTE:  The  update  rules  for  Ak  in  the  algorithm  can  be  found  in  (1)  and 
have  important  implications  for  the  convergence  of  the  algorithm. 

Figure  3.2  MGPS-RS  Algorithm  for  Stochastic  Optimization 
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The  following  assumptions  are  necessary  for  the  convergence  theory  of  the  MGPS-RS 
algorithm  to  hold  (50): 

A1  All  iterates  Xk  produced  by  the  MGPS-RS  algorithm  lie  in  a  compact  set. 

A2  The  objective  function  /  is  continuously  differentiable  with  respect  to  the  continuous 
variables  when  the  discrete  variables  are  fixed. 

A3  For  each  set  of  discrete  variables  Xd ,  the  corresponding  set  of  directions  Di  =  GiZi , 
as  defined  in  (2.1),  includes  tangent  cone  generators  for  every  point  in  ©c. 

A4  The  rule  for  selecting  directions  D\,  conforms  to  0C  for  some  e  >  0. 

A5  For  each  q  =  1,2 ,...,nc,  the  responses  {Fgs}^!=1  are  independent,  identically  and 
normally  distributed  random  variables  with  mean  f(Xq)  and  unknown  variance  aq  < 
oo,  where  aj  /  aq  whenever  £  ^  q. 

A6  The  sequence  of  significance  levels  { ar }  satisfies  ar  <  oo,  and  the  sequence  of 

indifference  zone  parameters  {<5r}  satisfies  lim,.—^  5r  =  0. 

A7  For  the  rth  R&S  procedure  considering  candidate  set  C  =  {Yi,  Y2,  ■  ■  ■ ,  Ync},  Pro¬ 
cedure  RS((7,  ar,  5r)  guarantees  correctly  selecting  the  best  candidate  Yrn  G  C 
with  probability  of  at  least  (1  —  ar)  whenever  /(Yfj)  —  /(Ym)  >  Sr  for  any  q  G 
{2,3, . . .  ,nc}- 

A8  For  all  but  a  finite  number  of  MGPS-RS  iterations  and  sub-iterations,  the  best  so¬ 
lution  Yni  G  C  is  unique;  i.e.,  f{Y mi)  /  /(Yrf/i)  for  all  q  G  {2,3,  ...,nc}  where 
C  =  {Yi ,  Y2, . . . ,  Yic}  C  M(X]t)  at  iteration  k. 

Sriver  includes  a  brief  discussion  on  these  assumptions.  Assumptions  Al,  A3,  and 
A5  are  noted  since  they  have  implications  in  the  model  formulation  and  presentation 
of  results  in  this  thesis.  Under  these  assumptions,  Sriver  et  al.  (52)  proved  almost 
sure  convergence  (i.e.,  with  probability  one)  of  a  subsequence  of  iterates  to  a  first-order 
stationary  point.  First-order  stationarity  in  a  mixed  variable  domain  was  first  formally 
defined  by  Lucidi  et  al.  (40)  in  the  context  of  a  more  general  framework  for  mixed  variable 
optimization. 
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Since  the  NOMADm  optimization  code  is  used,  it  is  important  to  highlight  the  spe¬ 
cific  R&S  procedure  implemented  in  the  code.  Based  on  the  results  of  the  computational 
results  in  (50),  Dunlap  (24)  chose  the  SSM  procedure  of  Pichitlamken  and  Nelson  (44) 
for  its  superior  performance.  SSM  is  a  fully  sequential  procedure  designed  for  use  with 
iterative  search  routines,  such  as  the  class  of  pattern  search  methods.  The  procedure 
collects  one  sample  at  a  time  from  every  candidate  and  eliminates  clearly  inferior  alter¬ 
natives  as  sampling  progresses.  To  identify  inferior  candidates,  SSM  performs  a  pairwise 
statistical  test  at  every  iteration  for  each  of  the  candidates.  By  removing  the  candidates 
that  are  statistically  inferior  to  all  the  other  members  of  the  candidate  set,  SSM  provides 
a  computationally  efficient  procedure  to  select  the  best  candidate.  Dunlap  states  that  an¬ 
other  advantage  of  SSM  is  the  ability  to  use  previously  stored  sampling  data  (24).  This 
mitigates  some  of  the  cost  of  obtaining  additional  simulation  responses  at  each  iteration 
of  the  optimization  algorithm. 

3.2  Multi-Echelon  Repair  Model  Formulation 

Since  the  multi-echelon  repair  problem  has  been  solved  analytically  under  certain 
assumptions,  these  solutions  can  be  used  as  a  baseline  to  compare  algorithm  performance  on 
similar  models.  The  specific  model  and  analytical  solution  used  as  a  baseline  in  this  thesis 
is  from  Gross  et  al.  (29),  who  gives  an  exact  solution  to  the  multi-echelon  problem  with 
exponential  distributions  for  holding  rates  and  arrivals  according  to  a  Poisson  process;  these 
results  are  used  as  a  baseline  comparison.  Gross  et  al.  devise  a  simple  two-echelon  model 
with  a  single  base  and  single  depot,  where  Al  machines  generate  item  failures  according  to 
a  Poisson  process  with  mean  A  =  pLjj.  A  failure  is  removed  from  the  machine  and  replaced 
with  a  spare  from  the  base  inventory.  If  no  spare  is  available,  the  machine  waits  for  an  item 
to  become  available,  either  through  repair  at  the  base  or  arrival  from  the  depot.  Failed 
items  extracted  from  a  machine  are  determined  to  be  either  base-serviceable  or  needing 
advanced  repair  that  cannot  be  performed  at  the  base.  A  certain  percentage  of  failed 
items,  a,  are  deemed  to  be  base  serviceable,  with  the  remaining  (1  —  a)  sent  to  the  depot 
for  repair. 
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Figure  3.3  Network  for  a  two-echelon  repairable  system 

Both  the  base  and  depot  repair  items  according  to  an  exponential  distribution  with 
mean  nB  and  y,Bl  respectively.  Each  base  has  a  specific  number  of  repair  channels.  An 
item  that  arrives  at  a  location  with  all  of  the  channels  busy  waits  in  a  queue.  There  is  a 
single  queue  at  each  location,  which  follows  a  FIFO  (first  in,  first  out)  rule.  When  a  repair 
is  complete,  the  item  is  returned  to  the  spares  inventory  unless  an  inoperable  machine  is 
awaiting  the  item,  in  which  case,  the  item  is  installed  and  the  machine  returns  to  service. 
However,  a  percentage  of  items,  (3 ,  after  undergoing  repair  at  the  base,  need  further  repair 
and  are  sent  to  the  depot. 

Recall  that  Gross  et  al  (29),  with  the  assumption  that  all  processes  are  Markovian, 
develop  an  expression  for  availability  in  terms  of  the  decision  variables.  The  decision 
variables  are  number  of  spares  y  to  be  inventoried  at  the  base,  and  the  number  of  channels 
at  the  base  and  depot,  cB  and  cB,  respectively.  The  number  of  spares  at  the  depot  is 
assumed  to  be  limitless.  A  diagram  of  the  flow  model  is  shown  in  Figure  3.3.  The  objective 
function  is  defined  to  be  a  linear  combination  of  the  costs  of  the  repair  channels  and  the 
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base-inventoried  spares.  Gross  et  al  (29)  define  the  problem  mathematically  as 

min  Z  =  kyy  +  kscs  +  koCD  (3.2) 

M+y 

S.t.  ^  pn  >  A, 
n=M 

where  M  is  the  number  of  machines  and  ki  is  the  cost  per  unit,  i  =  y,  B ,  D.  The  constraint 
is  absolute  availability  A,  with  pn  equal  to  the  steady  state  probability  that  n  units  are 
operational.  Gross  et  al  (29)  use 

M  =  5,ky  =  20,  ks  =  8,ko  =  10,  A  =  0.9,  a  =  0.5,  f3  =  0.5,  pjj  =  1,  fiB  =  llD  =  5 

as  the  parameters  for  this  model.  Unless  noted  otherwise,  these  parameter  values  are  used 
in  the  remainder  of  the  construction  of  the  simulation  and  mathematical  models. 

In  (3.2),  availability  is  the  stochastic  component  of  the  problem  to  be  supplied  by 
the  simulation  model.  However,  the  availability  level  appears  in  the  inequality  constraint 
and  not  in  the  objective  function.  Since  MGPS-RS  is  not  currently  able  to  handle  stochas¬ 
tic  nonlinear  constraints,  the  problem  formulation  given  in  (3.2)  must  be  modified.  Two 
approaches  for  doing  so  are  considered:  1)  reformulate  (3.2)  as  an  unconstrained  problem 
with  a  term  added  to  the  objective  function  to  penalize  infeasibility,  or  2)  swapping  the 
constraint  and  objective  function.  The  first  approach  becomes  an  unconstrained  mini¬ 
mization  problem  with  both  deterministic  and  stochastic  components.  In  the  latter  case, 
the  problem  becomes  one  of  maximizing  availability  (or  minimizing  expected  backorders, 
as  in  Sherbrooke  (48),  subject  to  a  linear  cost  or  budget  constraint.  This  problem  is  es¬ 
sentially  an  integer  knapsack  problem,  containing  a  linear  objective  function  subject  to  a 
single  linear  constraint  (56).  The  first  formulation  is 

minZ  =  20 y  +  8 eg  +  10cd  +  Pa{ 0.9  —  F),  (Formulation  A) 

where  Pa  is  the  penalty  and  F  is  the  absolute  availability  taken  from  the  simulation 
response.  Penalty  terms  are  well-suited  for  problems  in  which  some  of  the  constraint 
functions  require  noisy  response  evaluations  from  the  model,  since  it  cannot  be  determined 
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prior  to  simulation  if  a  design  is  feasible  with  respect  to  these  constraints.  However,  as 
in  the  deterministic  case,  it  is  noted  that  penalty  methods  may  suffer  from  computational 
difficulties  due  to  ill-conditioning  (50). 

The  second  formulation  is 


miwZ  =100 


1  - 


M+y  \ 

J2  Pn  ) 

,n=M  ) 


s.t.  20 y  +  8 cb  +  10cd  <  B 


(Formulation  B) 


where  B  is  the  budget  constraint.  The  availability  maximization  problem  is  converted 
to  a  minimization  problem  by  subtracting  from  1,  since  MGPS-RS  is  presented  as  a  min¬ 
imization  algorithm.  However,  the  problem  still  behaves  as  a  knapsack  with  an  increase 
in  the  decision  variables  yielding  a  better  result.  The  result  is  multiplied  by  the  constant 
100  simply  to  signify  a  percentage. 

To  apply  Formulation  A,  certain  parameters  must  be  determined.  The  first  is 
how  to  deal  with  values  above  the  availability  value  A.  In  this  formulation,  there  is  no 
cost  improvement  for  exceeding  the  availability  threshold.  Therefore,  the  penalty  is  only 
applied  for  values  of  availability  less  than  90%  and  all  responses  above  are  set  to  the  value 
of  0.9.  Also,  the  penalty  coefficient  for  the  new  objective  function  must  be  determined. 

Since  availability  was  originally  a  constraint,  it  is  important  that  the  penalty  coeffi¬ 
cient  in  the  objective  function  dominate  the  other  terms.  This  approach  attempts  to  force 
the  solution  to  satisfy  the  constraint,  and  is  consistent  with  the  general  penalty  function 
philosophy.  In  this  case,  the  coefficient  was  chosen  to  be  ten  times  larger  than  the  other 
terms.  Also,  the  square  root  of  the  deviation  from  the  constraint  is  taken  to  signify  that 
even  minor  violations  of  the  penalty  get  a  harsh  penalty.  Figure  3.4  shows  the  effect  of 
the  additional  penalty  term  on  the  objective  function.  The  resulting  model  is 

min  Z  =  20y  +  8 cB  +  10cD  +  2000^(0.9  -  F0,9)  (MODEL1A) 

y,  cB,  cD  e  z+ 
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I  1. 1  1.1  U  M  t.S  l.l  1.1  1.1  1.1 

Deviation  from  Availability  Constraint 

Figure  3.4  Penalty  assessed  for  violation  of  availability  constraint. 


where 


Fq.9  = 


F  if  F  <  0.9 
0.9  if  F  >  0.9. 


The  second  approach  only  requires  determination  of  the  budget  B.  A  budget  of  96 
monetary  units  was  chosen  to  maintain  a  comparison  with  the  analytical  solution.  The 
resulting  mathematical  model  is 


s.t.  20 y  +  8 cb  4-  10cd  <  96 
y,cB,cD  <E  Z+. 


(MODEL1B) 
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The  decision  variables  in  MODEL1A  and  MODEL1B  are  all  discrete.  This  causes 
the  decision  space  to  be  the  set  of  user-defined  neighbors  because,  with  no  continuous 
variables,  the  mesh  is  an  empty  set.  The  discrete  model  is  still  used  for  analysis  of  the 
algorithm  and  the  results  are  reported  in  Chapter  4.  To  test  a  true  mixed  variable  problem, 
the  model  is  further  adapted  to  provide  for  the  inclusion  of  continuous  variables.  In  the 
updated  model,  instead  of  the  number  of  repair  channels  being  decision  variables,  they  are 
fixed  at  three  each.  To  replace  them,  continuous  variables  representing  the  amount  of 
workday  the  base  and  depot  operate  are  introduced.  All  other  variables  and  parameters 
remain  the  same.  The  continuous  variables  can  be  likened  to  shift  work  or  overtime, 
for  example,  where  33%  indicates  a  location  is  open  eight  hours  each  day.  Changing 
the  model  to  a  mixed  variable  problem  eliminates  any  potential  for  comparison  to  the 
analytical  model  presented  in  (29). 

The  mixed  variable  problem  using  Formulation  A  is 


min  Z  =  20y  +  80.tb  +  100.td  +  2000^(0.9  -  F0,9)  (MODEL2A) 

xb,xd  €  [0.33, 1.00] 
y  £  Z+  and  xb,xd  £  M 


where 


Tq,  9  = 


F  if  F  <  0.9 
0.9  if  F  >  0.9. 


The  coefficients  for  the  base  and  depot  variables  increased  to  80  and  100,  respectively,  to 
accommodate  the  change  in  scale.  The  mixed  variable  problem  using  Formulation  B  is 
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(MODEL2B) 


s.t.  20 y  +  80 xb  +  IOOx’d  <  200 
xb,xd  G  [0.33, 1.00] 
y  G  Z+  and  xb,xd  G  K. 


Finally,  the  simulation  model  was  modified  to  incorporate  the  option  to  use  contract 
support  at  the  base.  Contract  support  is  a  fixed  cost,  but  if  used,  the  facility  can  repair 
items  quicker,  and  fewer  items  need  to  be  sent  to  the  depot  for  additional  repair  after 
attempting  repair  at  the  contractor  location.  The  validity  of  this  option  can  be  seen  in 
the  area  of  contractor  support  or  outsourcing.  The  use  of  contractor  support  is  repre¬ 
sented  by  a  binary  variable  w,  with  w  =  1  indicating  its  use.  The  simulation  model  was 
augmented  with  a  module  representing  contract  support  that  was  enabled  only  when  the 
binary  variable  was  set  to  one.  The  model  formulation  that  includes  the  contract  option 
is 


s.t.  20 y  +  80 xb  +  100xd  +  50 w  <  240 
xb,xd  G  [0.33, 1.00] 

V  G  Z+;  ib,idGM;  w  G  {0, 1}. 


(MODEL3) 


The  contracting  option  was  only  incorporated  into  the  budget  constrained  problem. 


3.3  Simulation  Model  Construction  and  Integration 

The  repair  model  was  constructed  using  the  widely  available  Arena1  discrete-event 
simulation  software.  Arena  is  easy  to  use  and  verify,  and  it  allows  input  variables  that 

xArenaw  simulation  software  is  a  registered  trademark  with  Rockwell  Software,  Inc. 
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can  be  specified  by  the  MGPS-RS  algorithm.  Also,  Arena  is  packaged  with  statistical 
tools  that  can  be  used  to  provide  the  necessary  availability  response  for  the  multi-echelon 
repairable  system. 

Steady-state  results  for  the  availability  are  desired.  However,  it  is  also  important  to 
have  a  degree  of  noise  in  the  availability  response  to  properly  test  the  algorithm.  With 
little  to  no  variability,  the  problem  becomes  too  much  like  the  deterministic  setting.  To 
balance  both  of  these  criteria,  each  simulation  replication  must  be  run  for  an  appropriate 
length.  The  approach  to  determine  the  simulation  length  is  similar  to  the  one  given  by 
Chrissis  and  Gecan  (19).  Ten  replications  were  divided  into  two  batches  of  five  each. 
The  five  replications  of  the  first  batch  were  run  for  a  specified  length  of  time  while  the 
other  five  are  run  for  a  slightly  longer  length.  Using  the  tools  provided  in  Arena’s  Output 
Analyzer,  a  hypothesis  test  to  compare  means  of  the  two  batches  was  performed.  If 
the  result  of  the  test  indicated  the  means  were  statistically  equal,  then  the  estimator  for 
availability  was  formed  using  the  replications.  If  the  hypothesis  test  result  showed  differing 
means,  the  replication  lengths  were  increased  and  the  process  was  repeated.  However,  this 
approach  became  too  expensive  computationally  to  perform  for  each  function  evaluation. 
The  method  was  still  used  to  determine  an  appropriate  length,  but  it  was  done  only  once 
to  determine  a  target  replication  length.  The  remaining  simulation  runs  were  conducted 
using  that  replication  length.  The  time  required  to  perform  one  replication  of  the  actual 
Arena  simulation  turned  out  to  be  negligible,  compared  to  the  cost  of  integrating  it  with 
the  optimization  code  and  the  setup  cost  involved. 

Integration  of  the  simulation  and  the  optimization  algorithm  was  particularly  diffi¬ 
cult  because  neither  the  optimization  code  nor  simulation  model  are  self-contained.  The 
NOMADrn  optimization  software  runs  inside  the  MATLAB  computing  environment  and 
makes  function  calls,  expecting  a  black  box  function  to  return  a  response.  However,  the 
simulation  model  runs  inside  the  Arena  simulation  software,  and  it  does  not  naturally  lend 
itself  to  easy  integration  with  a  MATLAB  code.  To  overcome  this  problem,  an  interface 
was  written  in  Microsoft  Visual  Basic  to  integrate  the  simulation  and  optimization  code. 
The  MATLAB  code  outputs  the  input  vector  and  waits  for  a  signal  from  Visual  Basic  to 
resume.  The  simulation  reads  the  input  vector,  runs  a  replication  given  those  parameters, 
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outputs  the  availability  response,  and  waits  for  the  next  input  vector.  The  simulation 
is  not  stopped  at  this  point  because  the  random  number  streams  would  reset,  meaning 
subsequent  identical  input  vectors  would  yield  the  exact  same  response,  resulting  in  no 
variability  or  noise  in  the  response. 

The  NOMADrn  software  website  (1)  gives  specifics  on  how  to  properly  define  the 
initial  iterate,  function  evaluation,  discrete  neighbors,  and  linear /bound  constraint  files  in 
the  software.  Before  the  algorithm  could  be  implemented,  certain  files  had  to  be  set  up 
which  define  the  problem  to  be  solved.  Among  them,  the  set  of  discrete  neighbors  was 
constructed  to  be  the  commonly  used  set  given  in  (2.4).  However,  as  described  in  Chapter 
2,  neighbors  for  categorical  variables  may  include  any  possible  value  for  the  categorical 
variable.  Since  the  contracting  variable  was  binary,  this  definition  of  the  neighbor  set  was 
sufficient.  The  linear  bound/linear  constraints  definition  was  important  for  two  reasons. 
The  first  was  to  ensure  assumption  A1  was  satisfied  that  the  iterates  lie  on  a  compact  set. 
Since  the  continuous  variables  x b  and  X£>  have  both  lower  and  upper  bounds,  A1  was 
satisfied.  The  second  reason  is  that  for  proper  definition  of  the  linear  constraints  is  that 
the  set  of  directions  must  positively  span  the  tangent  cone  of  the  nearby  linear  constraints. 

The  MGPS-RS  algorithm  within  the  NOMADrn  software  was  applied  to  the  multi¬ 
echelon  repair  models  explained  in  this  chapter.  The  results  for  each  of  the  models  are 
presented  in  the  next  chapter. 
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IV.  Results 


This  chapter  presents  the  results  for  each  of  the  models  outlined  in  Chapter  3.  For 
the  models  that  have  an  associated  deterministic  or  analytical  solution,  a  discussion  on 
accuracy  is  presented.  Also  given  is  a  brief  examination  regarding  the  effect  on  convergence 
of  the  choice  of  indifference  zone  and  alpha  parameters,  and  the  associated  decay  rates  used 
in  the  R&S  procedure. 

4-1  Discrete  Variable  Models 

In  the  case  with  all  discrete  variables,  the  POLL  step  of  the  MGPS  algorithm  simply 
consists  of  repeated  evaluations  of  discrete  neighbor  points,  and  the  extended  poll  step 
is  empty.  In  the  discrete  neighbor  search,  the  algorithm  acts  similar  to  a  one-factor-at-a- 
tirne  approach;  it  methodically  searches  the  neighbors  looking  at  each  of  the  variables  for 
the  best  improvement  by  increasing  each  subsequently  by  one. 

The  algorithm  was  run  using  both  models,  MODEL1A  and  MODEL1B,  at  two 
different  levels  of  machines,  M  =  5  and  M  =  20.  For  both  models  with  M  =  5,  complete 
enumeration  of  the  feasible  design  space  occurred.  This  is  primarily  due  to  a  small 
decision  space  and  each  of  the  variables  having  a  large  and  nearly  equal  effect  on  the 
objective  function.  However,  the  optimal  solution  was  obtained  in  both  cases  showing 
that  the  noise  in  the  simulation  response  was  overcome.  In  addition  to  the  analytical 
solution,  the  discrete  model  results  can  also  be  compared  to  the  proposed  particle  swarm 
optimization  approach  (PSO)  by  Alkhamis  and  Ahmed  (6).  It  should  be  noted  that  the 
PSO  algorithm  attempts  to  solve  the  original  formulation  (3.2)  and  Alkhamis  and  Ahmed 
do  not  provide  any  measure  of  computational  time.  Table  4.1  summarizes  the  results  for 
both  models  with  M  =  5  compared  with  the  analytical  solution  from  (29)  and  PSO. 

In  the  case  with  M  =  20,  the  results  for  the  discrete  models  differed.  The  budget 
constrained  availability  model,  MODEL1B,  behaved  similarly  to  the  problem  with  fewer 
machines.  However,  the  search  method  was  able  to  eliminate  some  inferior  solutions, 
so  the  solution  space  was  not  completely  enumerated.  Conversely,  MODEL1A  found 
no  improvement  from  the  initial  iterate  with  M  =  20  because  the  polling  of  the  discrete 
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Table  4.1  Model  1A/1B  Results,  M=5. 


y 

XB 

xD 

Evals 

True  (Analytical) 

3 

2 

2 

- 

MODEL1A 

3 

2 

2 

60 

MODEL1B 

3 

2 

2 

87 

PSO 

3 

2 

2 

n/a 

Table  4.2  Model  1A/1B  Results,  M=20. 


y 

XB 

XD 

Evals 

True  (Analytical) 

8 

5 

6 

- 

MODEL1A 

1 

1 

1 

16 

MODEL1B 

8 

5 

6 

760 

PSO 

9 

4 

6 

n/a 

neighbors  failed  to  yield  an  improved  solution.  This  can  be  attributed  to  the  cost  of 
an  increase  in  any  of  the  decision  variables  outweighing  each  of  their  contributions  to 
the  improvement  in  availability.  A  higher  penalty  on  the  deviation  from  the  desired  0.9 
availability  may  have  remedied  this  problem.  For  example,  an  initial  point  of  {2,  2,  2} 
when  M  =  5  would  have  likely  eliminated  reaching  {1, 1, 1},  provided  the  variability  in 
the  model  was  not  too  large.  Table  4.2  summarizes  the  results  for  the  discrete  variable 
models  with  M  =  20  compared  with  the  respective  analytical  solution  and  PSO. 

It  is  important  to  note  that  a  better  initial  design  vector  with  the  variables  at  higher 
starting  levels  would  have  likely  prevented  some  of  the  enumeration,  since  the  knapsack 
problem  with  negative  coefficients  yields  a  better  solution  relative  to  an  increase  in  the 
values  of  the  design  variables. 

4-2  Mixed-Variable  Models 

Results  for  MODEL2A  were  more  interesting.  The  size  of  the  problem  enabled  the 
continuous  variables  to  be  “discretized”  in  intervals  of  0.025,  creating  a  complete  discrete 
decision  space  with  a  finite  number  of  points.  A  deterministic  model  was  constructed  and 
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Table  4.3  Model  2A  Results. 


y 

XB 

xD 

Z 

MODEL2A 

5 

0.833 

0.9583 

412.494 

Deterministic 

2 

0.925 

0.95 

269 

run  using  all  possible  decision  points,  completely  enumerating  the  discrete  space,  to  give 
an  estimated  value  of  the  objective  function  without  the  random  noise.  The  algorithm 
was  unable  to  obtain  the  global  optimum.  The  results  for  MODEL2A  are  also  in  Table 
4.3;  the  deterministic  model  indicates  an  approximate  “optimal”  solution. 

A  problem  with  using  the  penalty  approach  is  that  /  is  not  continuously  differen¬ 
tiable,  violating  the  assumptions,  listed  in  Section  3.1,  necessary  to  ensure  convergence 
to  a  stationary  point.  Also,  the  algorithm  has  a  significant  problem  at  the  0.9  boundary 
point  particularly  because  of  the  random  noise.  Since  the  solution  is  known  to  have  a 
true  response  of  0.9041  (just  above  0.9)  from  the  deterministic  function,  noise  can  cause 
a  significant  penalty  to  be  assessed.  There  is  no  noise  in  the  objective  function  on  the 
lower  end  because,  if  the  availability  is  above  0.9,  the  objective  function  is  only  made  up 
of  deterministic  values.  Obviously,  the  error  in  the  objective  function  at  points  that  yield 
availability  values  close  to  0.9  is  not  normally  distributed.  In  fact,  it  is  only  one-sided. 
Figure  4.1  shows  the  behavior  of  the  objective  function  for  the  deterministic  model  in  terms 
of  the  continuous  variables  when  the  number  of  spares  is  two.  Notice  the  surface  sharply 
flattens  in  the  area  around  the  optimum,  indicating  the  availability  above  0.9.  Surface 
plots  for  different  number  of  spares  are  located  in  the  appendix. 

The  algorithm  had  no  problem  handling  MODEL2B.  Since  the  true  objective 
function  is  linear  with  a  single  linear  constraint,  the  objective  function  is  simply  a  series  of 
planes  that  are  “cut-off”  by  a  plane  representing  the  constraint.  An  approximate  solution 
was  found  in  less  than  750  function  evaluations.  An  estimated  true  optimal  solution  was 
also  obtained  for  this  model  using  the  discretized  deterministic  model.  Table  4.4  compares 
the  results  of  the  deterministic  model  with  MODEL2B. 
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Table  4.4  Model  2B  Results. 


y 

xB 

xD 

Z 

MODEL2B 

2 

0.8864 

0.8958 

13.95 

Deterministic 

2 

0.875 

0.9 

13.89 

Table  4.5  Model  3  Results. 


y 

w 

XB 

XD 

Z 

MODEL3 

2 

“Yes” 

0.3333 

0.8333 

25.69 

Results  were  similar  for  MODEL3  since  the  true  objective  function  in  MODEL3 
is  similar  to  MODEL2B.  The  algorithm  seemed  to  handle  the  binary  variable  quite  well. 
The  indicated  optimum  found  by  the  algorithm  is  shown  in  Table  4.5. 

4  ■  3  Indifference  Zone 

The  R&S  procedures  used  in  this  algorithm  require  the  user  to  explicitly  set  the 
initial  values  for  the  indifference  zone  5  and  significance  level  a,  and  the  respective  decay 
factors.  Sriver  (50)  does  not  give  guidance  for  the  specific  values  to  use,  only  that  the 
values  be  “loose.”  The  reason  5  must  be  initially  started  high  is  because  this  causes  more 
samples  to  be  required  in  the  beginning  of  the  algorithm  when  there  is  little  information 
known  regarding  the  variability  of  the  simulation  model.  As  additional  samples  are  taken, 
the  indifference  zone  may  be  reduced  to  identify  inferior  solutions.  The  results  of  changing 
the  indifference  zone  value  and  decay  rate  verify  this  assertion.  When  the  indifference  zone 
or  the  decay  rate  were  set  very  low,  the  algorithm  yielded  much  poorer  final  solutions  and 
required  many  more  function  evaluations.  The  default  values  used  by  Dunlap  (24)  in  the 
NOMADrn  software  performed  well;  however,  this  may  be  a  coincidence  with  the  scale 
of  the  problems  used  in  this  analysis.  For  MODEL2A,  the  indifference  zone  and  the 
indifference  zone  decay  were  each  varied.  The  results,  shown  in  Table  4.6,  show  that  a 
smaller  initial  indifference  zone  obtained  the  same  solution,  but  took  many  more  simulation 
runs.  The  smaller  decay  rate  yielded  worse  results  in  the  same  number  of  simulation  runs. 
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Table  4.6  Model  2 A  results  using  different  indifference  zone  parameters. 


MODEL2A 

Indifference  Zone 

IZ  Decay  Rate 

Z 

Evaluations 

Baseline 

100 

0.95 

412.49 

284 

Smaller  Indifference  Zone 

50 

0.95 

406.24 

500* 

Lower  IZ  Decay 

100 

0.75 

1028.87 

500* 

*Note:  The  number  of  function  evaluations  was  capped  at  500. 


Table  4.7  Model  2B  results  using  different  indifference  zone  values. 


MODEL2B 

Indifference  Zone 

Z 

Evaluations 

Baseline 

100 

25.73 

215 

Smaller  Indifference  Zone 

10 

26.53 

353 

A  similar  result  was  obtained  for  MODEL2B  using  a  smaller  indifference  zone. 
Table  4.7  shows  the  results. 

Finally,  it  should  be  noted  that  the  computational  limits  in  this  study  were  not  in  the 
optimization  code  nor  the  simulation  runs.  The  majority  of  the  overhead  came  from  the 
integration  and  transfer  of  data  between  the  optimizer  and  the  simulator.  The  advantages 
of  using  an  efficient  ranking  and  selection  procedure  are  usually  concerned  with  sampling 
cost,  the  number  of  simulations  runs  necessary  to  produce  a  “best”  solution  from  a  set  of 
candidates.  However,  sampling  cost  can  sometimes  be  outweighed  by  switching  cost,  or 
the  cost  of  switching  simulation  configurations.  Sriver  included  the  switching  cost  in  his 
analysis  of  R&S  procedures.  He  showed  that  switching  in  the  SSM  algorithm  can  be  quite 
significant,  “requiring  more  switches  than  SAS  by  approximately  two  orders  of  magnitude 
on  each  of  the  test  problems”  in  Sriver’s  computational  results  (50). 

The  results  presented  in  this  chapter  show  reasonable  results  for  the  use  of  the  MGPS- 
RS  algorithm  for  simulation  optimization.  The  next  chapter  provides  some  conclusions  of 
this  study,  as  well  as  recommendations  for  future  research  in  this  area. 
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V.  Conclusions  and  Recommendations 


This  chapter  summarizes  the  results  given  in  Chapter  4  and  gives  recommendations  for 
future  research.  Possible  future  work  includes  continued  analysis  of  the  algorithm  and 
research  on  integration  issues  between  simulation  and  optimization  applications. 

5. 1  Conclusions 

This  thesis  provided  a  framework  for  applying  the  MGPS-RS  algorithm  to  a  simula¬ 
tion  model  representing  a  real-world  system.  The  results  on  the  simple  two-echelon  model 
show  promise  that  the  algorithm  applies  well  to  simulation  optimization.  The  algorithm 
was  able  to  handle  random  noise  in  the  response,  moving  towards  improving  solutions. 
Additional  findings  show  the  selection  of  indifference  zone  parameters  may  have  a  signif¬ 
icant  effect  on  the  number  of  simulation  runs  needed  to  produce  a  solution  and  on  the 
quality  of  that  solution.  The  research  also  provided  possible  pitfalls  to  the  application  of 
the  algorithm  to  an  actual  simulation  environment.  Simulations  are  not  often  written  in 
the  software  that  contain  robust  optimization  algorithms,  so  the  need  for  software  integra¬ 
tion  becomes  critical,  particularly  when  the  number  of  simulation  runs  is  large.  Finally, 
the  NOMADrn  software,  the  only  software  that  implements  the  MGPS-RS  algorithm,  was 
effectively  integrated  with  an  external  simulation  application. 

5.2  Future  Algorithm  Research 

The  results  of  this  study  highlight  some  items  for  future  research.  The  first  is  to  de¬ 
velop  a  more  formal  design  and  analysis  than  used  in  this  thesis.  This  thesis  provides  only 
an  initial  study  regarding  the  use  of  this  algorithm  for  optimization  via  simulation.  Addi¬ 
tional  follow-on  computational  studies  should  have  a  formal  experimental  design  structure. 
Also,  a  study  on  a  larger  more  complex  model  may  yield  additional  insight  into  the  perfor¬ 
mance,  limitations,  and  abilities  of  the  algorithm.  The  results  from  the  tests  on  varying 
indifference  zone  parameters  warrants  an  examination  of  the  user-defined  R&S  parameters, 
to  include  the  choice  of  indifference  zone  and  its  decay  rate.  A  way  to  balance  sampling 
and  switching  cost  could  improve  the  performance  of  the  algorithm,  especially  when  inte- 
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gration  becomes  a  problem.  As  Sriver  suggested  (50),  a  simple  modification  may  be  to 
implement  the  minimum  switching  sequential  procedure,  developed  by  Hong  and  Nelson 
(31)  as  a  sequential  sampling  technique  that  uses  the  same  number  of  switches  as  two-stage 
procedures,  and  determine  if  its  performance  warrants  inclusion  in  the  algorithm.  With 
less  switching  of  the  input  parameters,  the  simulation  could  run  batches  of  replications 
before  returning  to  the  optimization  procedure. 

5.3  Future  Integration  Research 

Since  the  integration  of  the  optimization  code  and  the  simulation  model  provided 
the  greatest  cost  in  the  simulation  optimization  procedure,  future  research  should  focus 
on  mitigating  this  cost.  Sriver  discusses  switching  cost  in  his  analysis  of  three  proce¬ 
dures  considered  for  implementation  in  MGPS-RS  (50),  where  Rinott’s  procedure  incurs 
no  switching  cost,  SAS  calls  for  a  single  switch  for  each  candidate,  and  SSM  requires  a 
switch  each  time  an  additional  sample  is  needed  in  the  ranking  and  selection  procedure. 
Because  SSM  performed  well  with  GPS  methods  in  controlling  the  number  of  required 
function  evaluations  while  guaranteeing  the  proper  level  of  confidence  in  the  iterate  se¬ 
lection,  Dunlap  (24)  incorporated  SSM  into  NOMADrn.  For  the  software  to  be  effective 
in  situations  where  the  switching  cost  is  more  expensive  than  sampling  cost,  as  in  this 
thesis,  an  alternate  procedure  could  be  considered.  Finally,  consideration  should  be  given 
to  incorporating  the  optimization  code  and  the  simulation  model  into  the  same  software 
application.  For  example,  Boesel  and  Nelson  (15)  interface  their  genetic  algorithm  using 
ranking  and  selection  technique  with  AweSirn!  simulation  software  for  seamless  operation. 
Since  MGPS-RS  is  already  encoded  into  the  NOMADrn  MATLAB  software,  simulations 
using  Simulink/SimE vents  packages  from  MathWorks  might  be  considered.  However,  since 
these  MATLAB  packages  are  not  as  widely  available,  the  algorithm,  if  continued  testing 
shows  positive  results,  may  be  considered  to  replace  heuristics  in  more  popular  simulation 
software. 
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VI.  Appendix 

This  appendix  shows  surface  plots  for  the  objective  function  behavior  with  no  noise  for 

MODEL2A. 


Deterministic  Model2A,  Number  of  Spares  =  1 


Deterministic  Model2A,  Number  of  Spares  =  3 
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Deterministic  Model2A,  Number  of  Spares  =  4 


Deterministic  Model2A,  Number  of  Spares  =  5 
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